AD-A257  226 

lillllillilillllllllllliililli: 


Rensselaer 


!  •'  'n- 


October  22, 1992 


Dr.  Edwin  P.  Rood 
Scientific  Officer  Code:  1132F 
Office  of  Naval  Research 
800  North  Qioincy  Street 
Arlington,  VA  22217-5000 

Dear  Dr.  Rood: 


OTIC 

,*_SCTK 

(NQvasisee 

'  5 


The  purpose  of  this  letter  is  to  transmit  the  seventh  quarterly  report  for  ONR 
Grant  N00014-91-J-1271,  "An  Experimental  Study  of  Plunging  Liquid  Jet  Induced 
Air  Carryunder  and  Dispersion”  (Lahey  &  Drew  -  CoPI). 

This  report  period  has  been  primarily  concerned  with  the  development  of  a  two- 
fluid  Computational  Fluid  Dynamics  (CFD)  model  for  a  fi'ee  two-phase  jet.  This 
model  has  been  inplemented  into  PHOENICS  and  evaluated  on  a  RISC/6000  work 
station  and  the  CRAY  YMP-8  at  Stennis  Space  Center. 

We  first  present  a  state-of-the-art  two-fluid  model  for  air/water  flows,  which 
consists  of  the  phasic  three-dimensional  conservation  equations  of  mass  and 
momentum.  The  associated  closure  conditions  will  also  be  discussed,  along  wdth 
some  general  constraints,  followed  by  specific  constitutive  relations  for  low  void 
fraction  bubbly  two-phase  flows.  Finally,  the  resiilts  of  the  numerical  evaluation 
of  a  planar  plunging  liquid  jet  will  be  presented. 

For  isothermal,  incompressible  air/water  flows  the  appropriate  ensemble 
averaged  conservation  equations  for  a  free  two-phase  jet  are.  i  D]&i^WtrnON  STATEM 
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In  these  conservation  equations,  Ok  is  the  volume  fraction  corresponding  to 
phase-k,*  p  is  density,  y  is  velocity,  T  is  the  total  stress  tensor,  g  is  gravity  and  M^i 
is  the  interfacial  force  (per  unit  volume). 


Because  of  the  averaging  process  the  system  of  equations  have  more  unknowns 
than  equations.  In  order  to  have  a  closed  mathematical  system  one  must 
constitute  the  interfacial  momentum  transfers  and  the  Reynolds  stresses. 


One  simple  way  of  constituting  the  interfacial  momentum  transfers  is  by  For 
assuming  that  the  only  source  of  interfacial  momentum  interchange  is  drag.  The  i 
algebraic  drag  law  which  is  normally  used  is: 

Mi  =4CD^av 
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where  is  the  mean  bubble  diameter  and  Cp  is  the  drag  coefficient.  A  well 
known  drag  correlation  due  to  Wallis  is: 
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The  phasic  momentum  equations  of  the  resultant  two-fluid  model  are  given  by: 
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where  T^gj)  is  the  stress  tensor  corresponding  to  bubble-induced  pseudo¬ 
turbulence  and  ^si)  stress  tensor  corresponding  to  shear-induced 

turbulence.  The  later  stress  tensor  can  be  computed  using  the  k-e  model. 

A  clear  drawback  of  this  simple,  but  popular,  model  is  that  many  other 
mechanisms  of  momentum  interchange  are  not  taken  into  account.  For 
example,  the  so-called  virtual  mass  force,  the  lateral  lift  force,  the  interfacial 
pressure  distributions  and  the  turbulent  dispersion  terms  are  all  missing. 
Indeed,  it  is  now  well  known  that  it  is  necessary  to  use  a  model  that  takes  into 
account  these  mechanisms. 

A  more  complete  two-fluid  model  can  be  derived  using  "cell  model"  ensemble 
average  techniques.  This  model  is  rigorous  for  a  dilute  dispersion  of  spherical 
bubbles  in  an  inviscid  flow  with  weak  void  and  velocity  gradients.  The  resulting 
phasic  momentum  equations  are: 
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where,  the  virtual  mass  acceleration  is, 
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and,  for  spherical  noninteracting  bubbles: 
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In  order  to  numerically  evaluate  this  two-fluid  model  we  need  the  appropriate 
boundary  conditions.  There  is  no  general  theory  for  the  type  of  boimdary  and 
initial  conditions  that  a  system  of  partial  differential  equations  requires  in  order 
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to  have  a  unique  solution.  Nevertheless,  for  a  hyperbolic  system  Cauchy 
conditions  are  recommended. 

Most  researchers  in  the  past  have  evaluated  two-phase  models  as  initial  value 
problems  using  parabolic  numerical  techniques,  thus  it  was  not  necessary  to 
specify  boundary  conditions  at  the  outlet  of  the  integration  domain.  However, 
when  this  approach  is  used,  one  has  to  specify  the  pressure  distribution  in  the 
integration  domain.  For  many  cases  a  hydrostatic  pressure  distribution  is  a  good 
approximation.  Unfortunately,  when  a  parabolic  scheme  is  used,  one  cannot 
compute  flow  recirculation.  In  particular,  one  cannot  compute  gas  flow  reversal, 
an  important  feature  in  the  two-phase  jets  which  are  of  interest  in  this  study. 

In  this  work  we  have  used  an  elliptic  (ie,  boundary  value)  calculational  scheme. 
That  is,  we  have  numerically  evaluated  the  full  two-fluid  model  and  the 
associated  k-e  model  using  appropriate  boundary  conditions  at  the  inlet  and  the 
exit.  This  complicates  the  evaluation  procedure,  however  only  in  this  way  can  we 
obtain  a  prediction  of  gas  flow  reversal. 

Figure- 1  shows  the  integration  domain.  We  have  refined  the  grid  near  the 
symmetry  plane  of  the  planar  jet  because  the  gradients  are  the  steepest  there.  In 
the  axial  (i.e.,  z)  direction  we  have  increasingly  larger  cells  because  of  the 
decreasing  gradients.  The  boundary  conditions  which  have  been  used  are: 

INLET(z  =  0,-|<y<|) 


Gas  mass  flux  =  a  pg  V 

(7a) 

Liquid  mass  flux  =  (1  -  a)  V 

(7b) 
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It  shoiild  be  noted  that  the  boundary  conditions  associated  with  the  mass  balances 
are  only  required  at  the  inlet.  This  can  be  understood  based  on  the  fact  that  the 
mass  balance  is  a  first  order  partial  differential  equation.  Moreover,  given  the 
velocity  field,  the  equation  may  be  solved  for  the  void  fraction  along  the  model's 
characteristics.  With  the  void  fraction  given  at  the  inlet,  the  void  fraction  field  can 
he  readily  evaluated. 

Because  the  velocity  fields,  as  well  as  the  turbulent  kinetic  energy  and  the 
dissipation,  are  not  well  known  at  the  outlet,  we  have  specified  natural  boundary 
conditions.  That  is,  we  set  the  gradients  of  these  variables  in  the  flow  direction 
equal  to  zero.  In  our  case  the  outlet  boundary  conditions  were  found  to  have  a 
negligible  effect  except  for  the  last  two  rows  of  cells. 

The  numerical  results  presented  herein  correspond  to  a  planar  liquid  jet 
impacting  a  liquid  pool  using  the  two-fluid  model  given  by  Eqs.  (1)  and  (5).  The 
liquid  jet  velocity  at  the  location  of  impact  was  u^^  =  ^  ^  constant  bubble 

diauneter  of  2  mm  was  assumed.  The  initial  velocity  profile  was  uniform.  The  jet 
width,  h,  was  4.03  mm,  and  the  void  fraction  at  the  location  of  impact  was 
assumed  to  be  uniform  and  equal  to  5%.  The  inlet  turbulent  intensity  was  3%  and 
the  inlet  turbulent  dissipation  was  computed  using  Eq.  (7h).  The  inclination 
angle  of  the  jet  was  measured  with  respect  to  the  vertical  plane  (i.e.  9  =  0  degrees 
means  a  vertical  planar  jet).  The  integration  domain  had  an  extent  of  y  =  0.2  m  in 
the  lateral  direction  and  z  =  0.25  m  in  the  axial  direction.  The  k-e  model  for 
turbulence  used  the  constant  values  suggested  by  Launder  &  Spalding  for  single¬ 
phase  flow.  We  have  presented  results  corresponding  to  a  vertical  liquid  jet  (9  =  0 
degrees)  ujiless  otherwise  stated. 

Figure-2  shows  the  gas  and  liquid  axial  velocities  as  a  function  of  lateral  direction 
for  z  =  0.225  m  (note,  y  =  0.1  m  is  the  jet’s  plane  of  symmetry).  Both  velocity 
profiles  show  a  Gaussian-like  profile.  We  see  that,  due  to  buoyancy,  the  liquid 
velocity  is  always  higher  than  the  gas  velocity,  a  well-known  characteristic  of 


Dr.  Edwin  R.  Rood 
October  22, 1992 
Page  7 

downflows.  The  relative  velocity  was  approximately  0.3  m/s  which  is  very  close  to 
the  terminal  rise  velocity  of  the  bubbles. 

Figure-3  depicts  the  axial  liquid  velocity  as  a  function  of  lateral  position  for 
different  axial  positions.  The  curve  labeled,  z  =  0.00125  m,  is  right  under  the 
impacting  location  and  one  can  see  an  almost  imiform  velocity  profile  (u^  =  5  m/s). 
As  we  move  down  in  the  pool  the  jet  is  dispersed  due  to  momentum  interchange 
with  the  surrounding  fluid.  The  curve  z  =  0.225  m  is  the  same  one  shown  in 
Figure  2. 

Figure-4  shows  the  turbulent  kinetic  energy,  k,  as  a  function  of  lateral  position  for 
z  =  0.225  m.  The  curve  shows  the  characteristic  relative  minimum  in  k  at  the 
symmetry  plane. 

In  Figure-5  the  liquid  velocity  field  has  been  plotted  as  a  function  of  axial  and 
lateral  position.  The  length  of  the  arrows  is  proportional  to  the  liquid  velocity  at 
the  location  of  the  arrow's  tail.  The  arrow’s  tail  is  located  at  the  center  of  the 
computational  cell.  The  arrow  scale  is  in  the  lower  left  corner  (Uj,/  =  2  m/s).  The 

spreading  of  the  jet  can  be  easily  seen  in  this  plot  as  noted  previously.  Near  the 
location  of  jet  impact  (z  =  0)  the  axial  velocity  is  almost  uniform.  Because  of  the 
momentum  interchange  between  the  jet  and  the  surrounding  fluid,  liquid  is 
entrained  in  the  later^,  y,  direction.  Finally,  one  may  note  the  formation  of  two 
weak  recirculation  zones  near  the  y-boundaries  for  large  z. 

Figure-6  shows  a  contour  plot  of  the  axial  liquid  velocity.  The  lines  connect 
positions  with  the  same  axial  velocity  (equivelocity  lines).  The  outer  curve 
corresponds  to  u^^  =  0.25  m/s,  and  the  difference  between  successive  lines  is  0.25 

m/s. 

Figme-7  shows  a  contour  plot  of  the  void  fraction  for  one  half  of  the  jet.  The  outer 
line  connects  points  with  the  void  fraction  a  =  0.25%.  The  divergent  lines  show 
that,  when  only  drag  is  used  for  closure,  the  void  fraction  field  spreads  only 
slightly  as  z  increases. 

Figure-8  shows  a  vector  plot  of  the  liquid  velocity  field  for  an  inclined  planar  jet 
(0  =  45°).  We  rotate  the  integration  domain  45°  in  order  to  have  the  plane  y  =  0 
aligned  with  the  jet  orientation.  This  was  done  for  computer  time  economy 
purposes  and  in  order  to  minimize  numerical  diffusion.  It  can  be  seen  that,  as 
expected,  the  gas  drags  the  liquid  away  from  the  centerplane. 

Figure-9  shows  a  vector  plot  of  the  gas  velocity  field  for  an  inclined  planar  jet 
(9  =  45°).  Of  particular  interest  are  the  results  shown  in  the  upper  right  comer, 
the  gas  velocity  (weighted  with  the  void  fraction)  at  the  y-boundary.  This  shows  a 
reversal  of  the  gas  flow  rate. 
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The  calculations  presented  herein  correspond  to  a  simple  two-fluid  model.  The 
calculations  corresponding  to  the  full  two-fluid  model  (ie,  using  Eqs.  (7))  are 
currently  under  way.  These  predictions,  and  the  plunging  liquid  jet  data  using  a 
planar  nozzle,  are  expected  to  be  available  during  the  next  report  period. 

If  you  have  any  questions  concerning  this  report,  please  contact  me  or  Professor 
Drew  at  your  convenience  [Lahey:  (518)  276-8579;  Drew;  (518)  276-6903]. 

Sincerely  yours, 


Dr.  R.T.  Lahey,  Jr. 

The  Edward  E.  Hood,  Jr/Professor  of  Engineering 
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Administrative  Grants  Officer 
Director,  Naval  Research  Laboratory 
Defense  Technical  Information  Center 
D.A.  Drew  (RPI) 

F.  Bonetto  (RPI) 
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Profiles  for  z  =  0.225  m 
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Fig.  7  Contour  Plot  of  the  Void  Fraction  (Half  of  the  Jet)  for  the  Vertical  Jet 


Pig.  9  Vector  Plot  of  the  Gas  Velocity  for  the  Inclined  Jet  (0  =  45°) 
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